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I.  INTRODUCTION 


The  discipline  of  Computational  Fluid  Dynamics  (CFD)  has  progressed  to  the  point 
that  both  Euler  and  Navier-Stokes  flow  solvers  can  now  be  used  to  support  the 
aerodynamic  design  process  in  a  practical  manner  [1  through  7].  These  solvers  are 
generally  structured,  unstructured,  or  hybrid  (that  is,  not  “gridless”)  in  nature.  As  such, 
they  require  the  creation  of  a  solution  grid  on  which  to  compute  the  flowfield.  Hence,  the 
construction  of  an  appropriate  grid  is  the  key  to  obtaining  valid  results.  However,  this  is 
currently  an  enigmatic  art  beyond  the  proficiency  of  most  practicing  aerodynamicists. 

The  primary  difficulty  in  using  CFD  at  a  practical  level  is  the  uncertainty  of  creating 
an  initial  grid  that  will  produce  a  sufficiently  accurate  solution  for  the  problem  at  hand. 
Currently,  one  must  rely  on  specialized  expertise  and  experience  to  create  an  adequate 
initial  grid.  In  addition,  grid  resolution  studies,  perhaps  supplemented  by  Richardson’s 
extrapolation  technique,  must  be  conducted  to  ascertain  the  “goodness”  of  a  particular 
grid  [8  through  10].  All  this  must  be  done  before  further  calculations  can  be  performed 
with  an  assured  level  of  accuracy.  However,  the  time  required  to  perform  such  “check 
out”  computations  is  simply  not  available  in  a  practical,  fast-paced,  design  development 
or  problem  analysis  scenario.  In  such  an  environment,  quick  turn-around  is  of  paramount 
importance,  and  it  is  imperative  that  an  appropriate  grid  be  generated  the  first  time. 

Another  point  to  consider  is  that  the  practical  application  of  CFD  to  aerodynamic 
design  is  often  focused  only  on  producing  adequate  force  and  moment  coefficients  - 
products  that  are  used  in  flight  performance  analyses,  Six-Degree-of-Freedom  (6-DOF) 
flight  simulations,  guidance  and  control  system  design,  and  High  Level  Architecture 
(HLA)  simulations.  In  such  cases,  it  is  not  necessary  to  determine  detailed  properties  for 
the  entire  flowfield — only  the  pressure  and  shear  forces  acting  directly  on  the  body  of 
interest.  For  these  kinds  of  situations,  accuracy  really  only  needs  to  be  maintained  near 
the  body  region(s)  of  interest,  thus  relaxing  any  requirement  for  stringent  precision 
throughout  the  entire  solution  field. 

Further,  the  level  of  accuracy  required  exerts  a  strong  influence  on  both  the  amount 
of  effort  spent  on  grid  construction  and  the  amount  of  Central  Processing  Unit  (CPU) 
time  needed  for  the  solver  to  reach  convergence.  In  other  words,  if  engineering  accuracy 
is  satisfactory,  then  less-refined,  “faster-running”  grids  can  provide  the  desired 
information  in  a  more  timely  and  less  expensive  manner. 

In  an  attempt  to  assist  the  aerodynamic  designer  in  harnessing  the  powerful  tool  of 
CFD,  an  initial  effort  is  made  to  develop  some  practical  grid  construction  rules  for 
incompressible,  laminar  flow  over  a  flat  plate.  While  this  particular  flow  geometry  has 
limited  direct  application  value,  it  does  provide  a  starting  point  for  other  plate-like 
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geometries,  such  as  wall  boundaries,  airfoil  surfaces,  and  missile  fins.  In  addition,  the 
guidelines  developed  will  be  soundly  anchored  by  the  Blasius  solution  for  laminar  flow 
over  flat  plates.  This  well-established,  classical  analysis  provides  very  accurate  flow  data 
against  which  to  test  and  tweak  any  proposed  rules  of  thumb. 

II.  METHODOLOGY 

A.  Overall  Approach 

The  perspective  taken  in  this  study  is  that  of  an  aerodynamic  designer  seeking  a 
direct,  clear-cut  and  practical  means  of  creating  suitable  CFD  grids.  As  such,  no  effort  is 
made  to  utilize  Richardson’s  extrapolation  technique  per  se  in  practice.  While  a  similar 
approach  will  be  used  initially  to  develop  the  guidelines,  the  intent  (in  application)  is  to 
avoid  the  additional  calculations  required  by  this  method.  Rather,  the  idea  is  to  quantify 
the  expected  level  of  accuracy  apriori  for  simple,  incompressible,  flat  plate  flow  grids. 

The  approach  then  is  to  first  select  a  suitable  flow  solver  and  then  determine 
appropriate  run  parameters  and  boundary  conditions.  Then  the  appropriate  size  of  the 
computational  domain  will  be  decided  and  an  efficacious  grid  point  distribution  function 
established.  Afterwards,  the  required  initial  grid  point  (or  cell  size)  spacing  will  be 
ascertained,  and  a  correlation  of  error  with  initial  grid  point  spacing  will  be  constructed. 
Next  the  minimum  number  of  grid  points  (or  cells)  in  each  coordinate  direction  will  be 
investigated.  The  resulting  guidelines  will  then  be  tested  for  a  case  outside  the  bounds  of 
the  cases  used  for  their  development. 

B.  Selection  of  Flow  Solver 

In  selecting  a  flow  solver  for  the  conduct  of  this  study,  it  should  be  noted  that 
numerous  valid  possibilities  exist.  However,  it  should  also  be  noted  that  the  various 
solution  algorithms  and  implementations  available  are  likely  to  produce  results  with 
somewhat  different  accuracies  for  the  same  flow  situation.  Given  this  circumstance  and  a 
desire  for  the  results  to  be  generally  meaningful,  it  becomes  clear  that  the  selected  solver 
should  be  widely  representative  of  those  currently  available.  Further,  it  should  be  robust, 
computationally  efficient,  and  easily  obtainable;  particularly  since  these  criteria  are  likely 
to  be  key  to  the  aerodynamic  designer  when  selecting  a  practical  CFD  tool. 

Such  a  tool  should  also  be  finite-volume  based  since  such  methods: 

(1)  inherently  conserve  mass,  momentum,  and  energy  for  each  grid  cell  (versus  finite 
difference  methods  which  have  issues  maintaining  conservation  in  three-dimensional 
flows),  (2)  permit  the  use  of  arbitrary,  non-smooth  meshes  (and  are  thereby  more  general 
than  finite-difference  techniques),  (3)  adapt  more  naturally  to  zonal  boundaries  (and 
consequently  handle  more  complex  geometries  than  finite-difference  approaches), 

(4)  maintain  second  order  accuracy  at  grid  boundaries  (versus  first  order  for  finite 
difference  techniques),  and  (5)  can  handle  geometries  with  large  wall  curvature  and/or 
severely  non-orthogonal  grid  boundaries  (versus  finite-difference  methods  which  have 
issues  with  these  features)  [11  through  13]. 
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A  solver  meeting  these  criteria  is  the  Wind  code  [14  through  18]  developed  by 
the  NPARC  (National  Project  for  Applications-oriented  Research  in  CFD)  Alliance  [19]. 
The  NPARC  Alliance  is  a  partnership  between  the  NASA  Glenn  Research  Center  (GRC) 
and  the  Air  Force  Arnold  Engineering  Development  Center  (AEDC),  but  it  also 
incorporates  other  U.S.  Government,  industry,  and  academic  associates.  Wind  embodies 
the  merged  capabilities  of  the  NASTD  code  [20  -  21)  (the  primary  flow  solver  at  the 
former  McDonnell  Douglas  facility  in  St.  Louis,  now  The  Boeing  Company  at  St.  Louis), 
the  NPARC  code  [22]  (the  original  solver  of  the  NPARC  Alliance),  and  the  NXAIR  code 
[23  -  24]  (the  primary  solver  used  at  AEDC  for  store  separation  problems).  As  such,  it  is 
a  Government-owned,  three-dimensional,  general-purpose  flow  solver  for  the  time- 
dependent,  compressible,  Euler  and  Reynolds-Averaged  Navier-Stokes  (RANS) 
equations.  Consequently,  the  use  of  Wind  for  this  study  should  ensure  the  general  utility 
of  the  results. 

C.  Boundary  Conditions 

To  assess  the  appropriateness  of  various  boundary  conditions,  an  initial  domain 
size  had  to  be  specified  within  which  to  perform  computations.  Although  the  required 
domain  size  will  be  determined  later,  it  was  crucial  that  the  initial  work  be  done  within  a 
sufficiently  large  region  to  avoid  possible  errors  due  to  closeness  of  the  boundaries.  It 
was  found  in  Hirsch  [25]  that  a  distance  of  50  chords  or  larger  between  an  airfoil  and  a 
boundary  is  not  uncommon.  So,  a  two-dimensional,  rectangular  domain  was  constructed 
(using  Gridgen  [26])  that  extended  50  plate  lengths  above  the  plate  and  50  lengths  in  both 
streamwise  directions  from  the  center  of  the  plate. 

In  addition,  a  basic  flow  scenario  had  to  be  defined.  For  purposes  of 
commonality  with  expected  follow-on  work  with  turbulent  flow,  the  experimental 
configuration  cited  in  Reference  27  was  chosen.  This  configuration  consisted  of  a  large 
flat  plate  on  which  measurements  were  made  in  a  low  speed  wind  tunnel.  The  region  of 
interest  selected  for  the  current  work  was  the  first  foot  of  the  plate,  and  the  slowest 
measured  velocity  (58  feet-per-second)  was  identified  as  being  appropriate  for  laminar 
flow.  For  these  conditions,  the  boundary  layer  thickness  was  computed  and  used  to 
estimate  the  initial  grid  spacing.  It  was  presumed  that  1/10  of  the  thickness  0. 1  feet  from 
the  leading  edge  would  provide  reasonable  resolution.  However,  in  order  to  enable  the 
grid  to  fit  within  the  available  computer  memory,  the  initial  spacing  was  increased  by  a 
factor  of  10  to  2.5x10 3  L,  where  L  =  the  length  of  the  region  of  interest,  in  this  case  1 
foot.  While  such  coarse  resolution  was  inadequate  to  resolve  the  flow  physics,  it  was  felt 
to  be  sufficient  to  compare  the  relative  merits  of  various  boundary  conditions. 

A  “Viscous  Wall,”  no-slip  boundary  condition  was  specified  on  the  entire  flat 
plate  (including  the  part  downstream  of  the  region  of  interest),  and  an  “Inviscid  Wall,” 
slip  condition  was  denoted  on  the  boundary  immediately  upstream  of  the  plate.  These 
specifications  remained  the  same  throughout  the  duration  of  this  investigation.  However, 
various  types  of  definitions  were  explored  for  the  other  boundaries.  After  numerous 
computations,  it  was  discerned  that  using  “Freestream”  one-dimensional,  characteristic- 
based  conditions  on  the  upstream  and  upper  boundaries  with  “Outflow”  extrapolation  on 
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the  downstream  boundary  worked  best.  The  “goodness”  of  the  calculations  was  assessed 
by  comparing  velocity  contours  throughout  the  field,  velocity  profiles  near  the  end  of  the 
region  of  interest,  and  the  skin  friction  distribution  along  the  plate. 

D.  Run  Parameters 

In  terms  of  run  parameters,  the  previously  mentioned  calculations  revealed  the 
fastest,  yet  physically  realistic,  solutions  were  produced  by  the  “Scalar  Implicit”  mode. 
This  setting  was  exercised  to  invoke  the  scalar  implicit  form  of  the  approximately 
factored  implicit  scheme  [20]  in  all  coordinate  directions.  The  default  would  have  been 
to  apply  a  block  implicit  algorithm  in  the  “viscous”  streamwise  direction.  This  action 
resulted  in  the  viscous  contributions  being  confined  to  the  right-hand-side  terms  of  the 
approximately  factored  equation  -  the  viscous  terms  within  the  factorization  were 
ignored.  Otherwise,  the  default  settings  proved  satisfactory.  The  other  defaults: 

(1)  implemented  a  second-order,  upwind-biased,  Roe  algorithm  modified  for  stretched 
grids  (with  a  TVD  [Total  Variation  Diminishing]  minmod  limiter  using  a  limit  factor  of 
2)  for  the  inviscid  terms,  (2)  employed  second-order,  central  differencing  for  the  viscous 
terms,  and  (3)  and  stipulated  that  CFL=1.3. 

E.  Parallelism  and  Modes  of  Communication 

In  addition  to  the  previous  assessments,  an  evaluation  was  made  of  parallelism 
in  Wind  to  determine  an  optimum  number  of  CPUs  to  use.  The  goal  was  to  reduce  the 
wall  clock  time  while  maximizing  CPU  usage  for  a  grid  comprised  of  multiple  zones.  It 
was  found  that  for  an  Origin  2000  computational  platform,  Wind  1 .0  was  90  percent 
parallelizable  [28]  in  the  “Scalar  Implicit”  mode,  and  98  percent  parallelizable  in  the 
“Block  Implicit”  mode.  Ten  CPUs  appeared  to  offer  a  “speedup”  (CPU  time  /  wall  clock 
time)  of  five  for  scalar  implicit  calculations  and  almost  nine  for  block  implicit 
computations. 

While  these  results  indicated  the  use  of  10  grid  zones  and  10  CPUs,  they  did 
not  suggest  how  to  structure  the  grid.  So,  an  analysis  was  performed  to  learn  the  most 
effective  manner  to  do  so.  After  numerous  variations,  it  was  discovered  that  10  vertical 
slices,  each  with  1/10  of  the  total  grid  cells,  produced  the  best  results  in  the  shortest  time. 
An  attempt  to  use  10  horizontal  slices  simply  “blew  up,”  perhaps  due  to  the  absence  of 
physical  boundaries  in  the  layers  above  the  plate.  Other  structures  with  small  internal  and 
large  external  zones  were  tried  but  did  not  work  as  well  as  the  vertical  slices.  It  turned 
out  that  using  10  vertical  slices  produced  a  better-than-expected  speedup  of  seven  over  a 
single  zone  computation. 

To  reduce  the  wall  clock  time  even  further,  a  comparison  was  made  between 
the  Indirect  and  Direct  Parallel  Virtual  Memory  (PVM)  processor  communication  modes 
that  Wind  uses.  This  was  done  for  a  10-zone  grid  with  10  CPUs.  It  was  observed  that  the 
default  Indirect  PVM  mode  only  yielded  a  CPU  to  wall  clock  speedup  of  three — while  the 
Direct  PVM  mode  maintained  the  afore-mentioned  speedup  of  seven,  a  decrease  in  wall 
clock  time  of  60  percent.  Consequently,  the  Direct  PVM  mode  was  retained  for  all 
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subsequent  computations.  It  should  be  noted,  though,  that  according  to  the  Wind  User’s 
Guide  [29]  the  Direct  PVM  mode  can  only  be  exercised  on  single  multi-processor 
systems. 


At  times  during  this  effort,  the  Wind  code  was  updated.  When  it  was,  a 
comparison  was  made  with  the  previous  version  to  ensure  the  codes  produced  identical 
results.  Then  an  assessment  of  speedup  capabilities  was  made  and  the  fastest  method(s) 
selected  for  use.  For  Wind  2.0,  it  was  observed  that  only  84  percent  parallelism  was 
achievable.  This  was  attributed  to  code  fixes  that  prevented  the  skipping  of  entire  Fortran 
DO  loops.  However,  for  Wind  3.0,  the  newly  (at  the  time)  implemented  Message-Passing 
Interface  (MPI)  mode  of  communication  was  observed  to  be  three  times  faster  than  Wind 
2.0  with  Direct  PVM.  Wind  3.0  with  MPI  also  produced  an  additional  CPU-to-wall- 
clock  speedup  of  3.4  times  that  of  Wind  2.0  with  PVM. 

III.  DEVELOPMENT  OF  GUIDELINES 


A  summary  of  the  boundary  conditions,  run  parameters,  and  processor 
communication  modes  used  to  develop  the  guidelines  is  presented  in  Table  1.  All 
associated  calculations  were  made  with  grids  composed  of  10  vertical  slices  using  10 
CPUs  as  described  above. 

Table  1.  Boundary  Conditions,  Run  Parameters,  and  Communication  Modes  Used  to 

Develop  Guidelines 


Boundary  Conditions 

Location 

Value 

Flat  Plate 

Viscous  Wall 

Stagnation  Streamline 

Inviscid  Wall 

Upstream  Boundary 

Freestream 

Upper  Boundary 

Freestream 

Downstream  Boundary 

Outflow 

Run  Parameters 

Parameter 

Value 

Solution  Methodology 

Scalar  Implicit 

Inviscid  Algorithm 
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Wind  1.0 

PVM  Direct 
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PVM  Direct 
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A.  Assessment  of  Domain  Size 


The  assessment  of  domain  size  began  by  using  conformal  mapping  to 
analytically  estimate  the  appropriate  location  for  freestream  boundaries  [30].  This 
approach  succeeded  in  producing  a  relationship  between  influence  distance  and  the 
accompanying  variation  of  velocity  from  freestream  conditions.  For  a  1  percent  variation 
in  freestream  velocity,  it  was  estimated  that  the  grid  boundaries  could  be  located  at  a 
distance  of  2.5L  from  the  center  of  the  plate. 

In  concert  with  this,  a  computational  appraisal  was  conducted.  This  evaluation 
began  with  the  “50L”  grid  described  in  the  previous  section.  It  continued  by  cutting  the 
distance  from  the  middle  of  the  region  of  interest  to  the  upstream,  downstream,  and  upper 
grid  boundaries  by  approximately  half  (to  the  closest  existing  grid  line)  each  time.  This 
was  done  to  avoid  shifting  the  internal  grid  lines  and  thereby  make  comparisons  with 
otherwise  identical  grids.  In  the  end,  grids  denoted  as  20.7L,  8.78L,  3.98L,  1.89L,  1.0L, 
0.73L,  0.54L,  and  0.50L  (to  designate  the  distance  to  the  boundaries  in  terms  of  the 
length,  L,  of  the  region  of  interest)  were  examined.  The  results  showed  that  the  boundary 
distance  for  this  incompressible  case  could  be  decreased  to  1.89L  with  no  detrimental 
effect  on  any  of  the  flow  variables  near  the  body.  If,  however,  it  had  been  necessary  to 
match  the  normal  velocity  component  far  from  the  plate,  the  boundary  distance  could 
only  have  been  decreased  to  20.7L.  In  fact,  the  normal  velocity  was  found  to  be  the  most 
sensitive  indicator  of  parameter  effects  throughout  the  entire  investigation. 

To  explore  the  previous  finding  further,  the  streamwise  (x-coordinate  direction) 
non-dimensional  velocity  component  of  the  50L  grid  was  examined.  It  was  scrutinized 
along  longitudinal  grid  lines  between  the  upstream  boundary  and  the  leading  edge  of  the 
plate  for  y-coordinate  locations  of  10'6L,  10'5L,  10'4L,  10'3L,  10'2L,  0.1L,  L,  and  10L 
units  above  the  plate.  It  was  noticed  that  the  velocity  component  remained  within 
0.25  percent  of  the  freestream  value  for  y  <  0. 1L  from  the  upstream  boundary  to 
x  <  -1.5L,  that  is,  one  plate  length  upstream  of  the  leading  edge.  This  was  also  true  for 
y  =  10L.  For  y  =  L,  however,  the  variation  was  larger,  approximately  0.5  percent,  at 
x  =  -6L.  While  this  variation  was  larger  and  extended  much  further  upstream,  it  was  still 
less  than  the  1  percent  variation  used  to  analytically  estimate  the  boundary  placement. 

From  these  results,  it  was  concluded  that:  (1)  the  freestream  boundaries  for  this 
incompressible  problem  could  be  located  2L  units  away  from  the  center  of  the  plate  with 
no  adverse  effect,  and  (2)  the  analytical  estimate  provided  a  valid,  yet  conservative 
estimate  of  freestream  boundary  placement. 

To  reduce  the  domain  size  even  further,  the  streamwise  boundaries  were  kept  at 
1.89L  while  the  upper  boundary  was  located  1.89L,  1L,  0.5L,  0.25L,  and  0.05L  units 
above  the  plate.  The  non-dimensional  streamwise  velocity,  u/Uoo,  was  unaffected  in  all 
cases,  yet  the  non-dimensional  normal  velocity,  v/Uoo,  was  affected  by  the  0.5L  and 
“shorter”  upper  boundary  placements.  This  consequence  explains  the  behavior  of  the 
Wind  validation  calculations  performed  by  Slater  [31].  His  results  agreed  well  with  the 
Blasius  solution  for  u/Uoo,  but  were  a  bit  low  for  v/Uoo.  Examination  of  his  upper 
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boundary  location  revealed  it  to  be  placed  at  y  =  0.08797L — well  below  the  distance 
found  to  influence  v/Uoo. 

In  any  case,  these  results  made  it  clear  that  a  1 .89L  (upstream)  by  1L  (upper)  by 
1.89L  (downstream)  domain  size  would  prove  adequate  to  analyze  incompressible  flat 
plate  flow.  However,  in  order  to  simplify  things  a  bit,  a  domain  size  of  2.0L  by  1.0L 
by  2.0L  was  used  for  the  following  analyses.  A  sketch  of  this  grid  and  the  applied 
boundary  conditions  is  shown  in  Figure  1 . 
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Figure  1.  2.0L  by  1.0L  by  2.0L  Computational  Domain 


B.  Grid  Point  Distribution  Function 

The  most  dominant  attribute  affecting  accuracy  is  the  spacing  of  grid  points  in 
the  solution  grid.  It  is  for  this  reason  that  analyses  and  studies  have  been  conducted  to 
discover  the  relative  accuracy  of  various  grid  point  distribution  techniques.  For  example, 
Vinokaur  [32]  performed  a  generalized  truncation  error  analysis  and  found  that  the 
hyperbolic  sine  (sinh)  function  produced  the  least  error  for  one-sided  point  clustering,  as 
would  be  typical  of  a  wall  to  freestream  grid.  He  further  found  that  the  hyperbolic 
tangent  (tanh)  function  produced  the  least  error  for  two-sided  clustering,  as  would  be 
required  for  a  wall-to-wall  grid. 
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Thompson  and  Mastin  [33]  quantitatively  examined  the  truncation  error  of 
10  specific  grid  point  distribution  functions  for  scenarios  of:  (1)  a  fixed  distribution 
function  with  a  varying  number  of  grid  points,  and  (2)  a  fixed  number  of  points  with 
varying  minimum  point  spacing.  In  general,  they  found  that  for  a  central  difference 
discretization,  the  solution  retained  second-order  accuracy  regardless  of  the  form  of  the 
point  distribution  function.  However,  only  4  of  the  10  distribution  functions  (the 
exponential,  hyperbolic  sine,  hyperbolic  tangent,  and  error  functions)  were  seen  to 
produce  suitable  grid  point  spacing  when  the  minimum  spacing  was  small,  as  would  be 
typical  near  a  viscous  wall.  Their  analysis  further  revealed  the  hyperbolic  sine  (sinh) 
produced  the  least  truncation  error  (of  these  four  functions)  near  the  wall,  while  away 
from  the  wall  the  hyperbolic  tangent  (tanh)  function  produced  less  error  than  either  the 
exponential  (exp)  or  hyperbolic  sine  functions.  The  error  function  (erf)  was  shown  to 
behave  similarly  to  the  hyperbolic  tangent  function  while  contributing  less  error  away 
from  the  wall.  However,  it  was  also  shown  to  add  substantially  more  error  near  the  wall. 
To  compare  the  two  most  promising  distribution  functions,  the  ratio  of  sinh  error  to  tanh 
error  was  examined  [30].  This  comparison,  made  throughout  the  entire  field  with  small 
values  of  minimum  point  spacing,  showed  that  the  hyperbolic  tangent  function  provides 
the  most  accurate,  general  purpose  distribution  of  grid  points. 

C.  Error  as  a  Function  of  Initial  Grid  Point  Spacing 

To  perform  this  phase  of  the  inquiry,  the  domain  shown  in  Figure  1  was 
“gridded  up”  using  an  initial  spacing  of  0.00 125L  (half  the  previous  value)  in  the  upper 
and  in  both  streamwise  directions,  centered  at  the  leading  edge.  It  was  acknowledged  that 
this  was  still  too  coarse  to  properly  resolve  the  flow  physics,  but  it  would  prove  useful  in 
comparing  the  relative  merits  of  various  grid  cell  distributions. 

Computations  were  first  made  for  a  uniform  grid.  After  that,  they  were  also 
made  for  grids  maintaining  the  same  uniform  streamwise  spacing  but  exercising  the  tanh 
point  distribution  function  in  the  y-direction.  These  grids  were  constructed  using  the 
same  initial  point  spacing  with  Vi,  !4,  1/10,  and  1/20  of  the  number  of  “uniform  grid” 
points  in  that  direction.  It  was  revealed  that  the  number  of  points  in  the  y-direction  could 
be  reduced  to  5  percent  of  that  required  for  a  uniform  grid  with  no  adverse  effect. 
Subsequent  use  of  the  tanh  function  in  both  streamwise  directions  (centered  at  the  leading 
edge)  proved  that  a  reduction  to  50  percent  of  that  required  for  a  uniform  grid  also  had  no 
adverse  effect. 

At  this  point,  the  initial  grid  spacing  was  halved  three  successive  times  and  the 
error  in  drag  coefficient,  Cd,  calculated  at  0.9308L  (a  point  where  measurements  were 
made  in  Reference  27.  In  addition,  the  number  of  streamwise  points  was  reduced  to 
25  percent  and  10  percent  of  the  uniform  value.  The  10  percent  case  was  found  to  have 
only  a  slight  effect  on  the  v/Uoo  profile,  but  not  Cd.  At  this  point,  it  is  worth  noting  that 
use  of  the  tanh  distribution  function  reduced  the  required  number  of  grid  points 
substantially  as  compared  to  a  uniform  grid.  The  reduction  to  10  percent  and  5  percent  in 
the  streamwise  and  normal  directions,  respectively,  diminished  the  total  number  to 
0.5  percent  of  uniform  number  while  maintaining  solution  accuracy. 


After  the  previous  computations  were  completed,  the  error  in  Cd  (as  compared 
to  the  Blasius  value)  was  correlated  with  initial  grid  spacing  [30].  The  resulting  function 
was  used  to  predict  the  initial  spacing  needed  to  keep  the  error  around  10  percent — an 
acceptable  level  of  accuracy  for  many  initial  and  intermediate  design  applications. 

The  predicted  spacing  was  3.36xl0"4L  units  which  was  applied  to  the  domain 
shown  in  Figure  1.  The  ensuing  grid  had  121  streamwise  points  upstream  of  the  plate, 
201  points  along  the  entire  length  of  the  plate  (171  of  which  lay  in  the  region  of  interest, 
L),  and  41  points  perpendicular  to  the  plate.  This  grid  produced  an  error  in  CD  of  only 
8.85  percent,  even  less  than  the  “requested”  value  of  10  percent.  This  performance 
demonstrated  the  guidelines  to  be  capable  of  providing  better-than-specified  accuracy. 
Figure  2  exhibits  the  convergence  curve,  that  is,  coefficient  versus  cycle  number,  for  Cd- 
Note  that  a  cycle  represents  a  complete  iteration  of  the  entire  grid.  Wind  computed  five 
iterations  in  each  zone  and  then  updated  all  zonal  boundary  data  to  complete  a  single 
cycle. 


Figure  2.  Cd  Versus  Cycle  for  Initial  Test  of  Guidelines 

While  the  guidelines  performed  better  than  expected  for  their  intended  purpose, 
Figures  3  through  5  show  they  functioned  similarly  well  for  the  velocity  profiles  and  skin 
friction  coefficient.  Figures  3  and  4  present  the  velocity  profiles  near  the  end  of  the 
region  of  interest  with  error  bounds  of  ±5  percent  and  ±10  percent.  It  is  clearly  seen  that 
the  computed  profile  for  streamwise  velocity,  u/U»,  overlays  the  Blasius  solution  while 
the  normal  velocity,  v/Uoo,  is  well  within  -5  percent  of  the  Blasius  curve.  Figure  5 
exhibits  the  computed  local  skin  friction  coefficient  overlaying  all  but  the  first  tenth  of 
the  Blasius  values.  This  difference  is  to  be  expected  since  the  grid  resolution  was  not 
intended  to  capture  the  physics  in  this  part  of  the  field  and  is  therefore  too  coarse  to  do  so. 
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v/Uinf 


Figure  3.  u/Ucc  Profile  at  x=0.9308Lfor  Initial  Test  of  Guidelines 
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Figure  4.  v/Um  Profile  at  x=0.9308Lfor  Initial  Test  of  Guidelines 
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Figure  5.  Local  Skin  Friction  Coefficient  Along  Plate  for  Initial  Test  of  Guidelines 

D.  Number  of  Grid  Points  in  Each  Direction 

With  the  legitimacy  of  the  guidelines  established,  an  endeavor  was  made  to 
further  reduce  and  quantify  minimum  numbers  of  grid  points  needed  in  each  coordinate 
direction.  The  previous  grid  was  evaluated  with  41,  31,  21,  16,  and  1 1  grid  points  in  the 
y-direction.  It  was  discovered  that  at  least  21  normal  points  were  required  to  stay  within 
an  accuracy  of  10  percent.  This  indicated  the  minimum  required  number  of  “y-direction” 
grid  points  to  be  approximately  1  percent  of  the  number  required  for  a  uniform  grid  [30]. 

In  the  streamwise  direction,  the  previous  grid  with  31  perpendicular  and  121 
upstream  points  was  used  with  variations  of  201,  151,  101,  51,  and  26  points  along  the 
plate.  This  produced  grids  with,  respectively,  171,  130,  88,  46,  and  24  streamwise  points 
in  the  region  of  interest.  The  downstream  grid  lines  were  allowed  to  vary  as  required  by 
the  tanh  distribution  function;  they  were  not  constrained  at  the  end  of  the  region  of 
interest.  This  examination  found  that  at  least  5 1  streamwise  points  (46  in  the  region  of 
interest)  were  needed  to  keep  the  error  within  10  percent.  This  yielded  a  minimum 
requirement  of  around  2  percent  of  the  number  required  to  fill  a  uniform  grid  between  the 
leading  edge  of  the  plate  and  the  downstream  boundary  [30]. 

Continuing  in  this  vein,  the  upstream  region  was  explored  using  the  previous 
grid  with  151  points  above  the  plate  and  31  normal  to  the  plate.  Variations  of  121,  91, 
61,31,  and  16  upstream  longitudinal  points  were  inspected  with  little  change  occurring  in 
the  error  for  Cq.  However,  “kinks”  occurred  in  the  velocity  profiles  when  less  than  3 1 
upstream  points  were  used.  Hence,  a  minimum  requirement  was  determined  to  be  about 
2  percent  of  the  number  required  to  uniformly  fill  the  region  between  the  upstream 
boundary  and  the  leading  edge  of  the  plate  [30]. 
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IV.  TESTING  THE  GUIDELINES 


To  this  point,  the  guidelines  had  only  been  tested  for  a  case  similar  to  those  from 
which  they  were  developed.  While  they  performed  quite  well,  they  needed  to  be  tested 
against  a  different  laminar  flat  plate  scenario  to  evaluate  their  generality.  The  instance 
selected  was  the  high  Reynolds  number  case  of  Schubauer  and  Klebanoff  [34]  that  is 
presented  in  Reference  35.  Although  the  flow  velocity  is  only  slightly  higher,  80  feet-per- 
second,  the  measurement  station  of  interest  is  located  at  5.25  feet  so  the  Reynolds  number 
is  an  order  of  magnitude  larger  at  2.8x1 06. 

The  grid  was  constructed  with  a  6-foot  region  of  interest  using  all  the  previously 
developed  guidelines.  Thus,  the  upstream  and  downstream  boundaries  were  placed 
12  feet  from  the  center  of  this  region  with  the  upper  boundary  6  feet  above  the  plate.  The 
initial  grid  spacing  was  predicted  to  be  6.9xl0'4  feet  (for  10  percent  accuracy)  and  the 
tanh  function  was  used  to  distribute  the  points.  Two  hundred  sixty  cells  (261  points) 
were  placed  in  the  upstream  region,  440  cells  (44 1  points)  above  the  plate,  and  87  cells 
(88  points)  perpendicular  to  the  plate.  The  results  proved  to  be  even  better  than  expected 
with  an  accuracy  in  Cd  of  9.6  percent  at  12,000  cycles,  2.6  percent  at  24,000  cycles,  and 
1.9  percent  at  36,000  cycles.  The  corresponding  convergence  curve  is  shown  in  Figure  6. 

— — —  Drag 


Cycle 

Figure  6.  Cd  Versus  Cycle  for  Reference  35  Case 


The  results  were  equally  good  for  the  velocity  profiles  and  local  skin  friction 
coefficient  which  are  presented  in  Figures  7  through  10.  Figure  7  illustrates  the  fact  that 
even  though  the  streamwise  profile  at  12,000  cycles  does  not  overlay  the  Blasius  solution 
(as  the  measurements  and  the  solutions  at  24,000  and  36,000  cycles  do),  the  computed 
drag  coefficient  still  has  an  accuracy  better  than  the  10  percent  value  that  was  specified. 
Figure  8  shows  that  the  solution  at  36,000  cycles  precisely  overlays  the  Blasius  curve  and 
does  not  observably  depart  from  it.  Figure  9  exhibits  the  normal  velocity  at  36,000  cycles 
as  nearly  overlaying  the  Blasius  curve  and  only  intaiding  slightly  into  the  +5  percent 
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region.  Lastly,  Figure  10  demonstrates  the  calculated  local  skin  friction  curves  at  24,000 
and  36,000  cycles  also  overlay  the  Blasius  values  except  near  the  leading  edge,  as 
expected.  Again,  this  figure  underscores  the  point  that  precise  capturing  of  the  field 
properties  is  not  necessary  to  obtain  a  useful  force  coefficient  for  design  purposes. 

-  12000  cycles 

—  —  24000  cycles 

-  36000  cycles 

-  - .  Blasius  Solution 
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Figure  7.  u/U «,  Profile  at  x= 5. 2 5ft  for  Reference  35  Case 
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Figure  8.  u/Um  Profile  at  x~ 5. 25ft  for  Reference  35  Case  After  36000  Cycles 
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Figure  9.  v/Ux  Profile  at  x=5. 25ft  for  Reference  35  Case  after  36000  Cycles 


-  12000  cycles 
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-  36000  cycles 
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Figure  10.  Local  Skin  Friction  Coefficient  Along  Plate  for  Reference  35  Case 


V.  SUMMARY 


The  ability  to  predict  or  specify  the  level  of  accuracy  for  a  CFD  grid  has  been 
demonstrated  for  incompressible,  laminar,  flat  plate  flows.  This  is  a  capability  that  is  not 
known  to  the  author  to  have  previously  existed.  In  addition,  it  has  been  proven  that  the 
required  size  of  the  domain  for  incompressible  flow  can  be  predicted  and  is  much  smaller 
than  would  otherwise  be  expected.  Further,  the  number  of  grid  points  required  to  obtain  a 
valid  solution  is  likewise  predictable  and  substantially  less  than  traditionally  thought 
provided  the  hyperbolic  tangent  function  is  used  to  distribute  them. 

Although  the  grid  generation  guidelines  that  provide  this  capability  were  formulated 
for  a  simple  type  of  flow,  they  provide  a  proven  starting  point  for  further  development. 
They  can  also  be  used  to  provide  insight  into  grid  construction  for  plate-like  geometries 
(such  as  flat  wall  boundaries,  airfoil  surfaces,  and  missile  fins)  in  low  speed  flow. 
Moreover,  they  can  be  used  directly  for  laminar  flow  problems  of  interest,  such  as  flows 
associated  with  biological  Micro-Electro-Mechanical  Systems  (MEMS)  devices.  In  short 
they  will,  as  intended,  assist  the  aerodynamic  or  fluid  dynamic  designer  in  harnessing  the 
powerful  tool  of  CFD. 
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